# -*- coding: UTF-8 -*-
import numpy as np
import matplotlib.pyplot as plt
# Pa 轴功率
Q = np.array([0.0,199.6,322.4,421.3,522.4,573.3,622.8,672.0,716.7,733.9,791.6,890.9,980.3])
Pa = np.array([133.86,185.33,218.08,249.18,282.56,298.04,315.25,331.08,344.46,349.81,367.13,396.67,423.18])
#yita = np.array( [0.0,45.4,62.3,71.3,76.9,78.7,79.6,80.3,81.1,81.1,81.5,80.5,78.1])

def gen_coefficient_matrix(X,Y):
    N = len(X)
    m =3
    A = []

    for i in range(m):
        a = []
        for j in range(m):
            a.append(sum(X**(i+j)))
        A.append(a)
    return A

def gen_right_vector(X,Y):
    N =len(X)
    m = 3
    b = []
    for i in range(m):
        b.append(sum(X**i * Y))
    return b
A = gen_coefficient_matrix(Q,Pa)
b = gen_right_vector(Q,Pa)

a0,a1,a2 = np.linalg.solve(A,b)

_X = np.arange(0,1000)
_Y = np.array([a0 +a1*x + a2*x**2 for x in _X])

plt.plot(Q,Pa,'ro', _X, _Y,'b',linewidth=2)
plt.ylim(ymax=500,ymin=0)
plt.title("y = {} + {}x + {}$x^2$".format(a0,a1,a2))
plt.annotate('Pa-Q', xy=(716.7,344.46), xytext=(800, 60), arrowprops=dict(facecolor='black', shrink=0.05))
plt.show()